library(tidyr)
library(plyr); library(dplyr)
library(foreign)
library(rdrobust)
library(arm)
library(lfe)
library(ggplot2)
library(xtable)
library(rdlocrand)

lag.new <- function(x, n = 1L, along_with){
 index <- match(along_with - n, along_with, incomparable = NA)
 out <- x[index]
 attributes(out) <- attributes(x)
 out
}

lead.new <- function(x, n = 1L, along_with){
 index <- match(along_with + n, along_with, incomparable = NA)
 out <- x[index]
 attributes(out) <- attributes(x)
 out
}


## Graphs
scale1sd <- function (x, center=TRUE) {
  s <- sd(x, na.rm=TRUE)
  m <- if (center) mean(x, na.rm=TRUE) else 0
  (x - m) / (1 * s)
}
scale2sd <- function (x, center=TRUE) {a
  s <- sd(x, na.rm=TRUE)
  m <- if (center) mean(x, na.rm=TRUE) else 0
  (x - m) / (2 * s)
}
mysw <- function (...) setwd(paste(..., sep = "/"))
mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  print(pdfname <- paste(paste(name, collapse=""), date, ".pdf", sep=""))
  pdf(pdfname, ...)
}
clusterSE <- function(x, method="arellano", cluster="group") {
  if (any(grepl("plm", class(x)))) {
    coeftest(x, plm::vcovHC(x, method=method, cluster=cluster))
  }
}
tri <- function (x, h, c=0) pmax(0, 1 - abs((x - c) / h))
is.even <- function(x) x %% 2 == 0
council_all<-read.dta("council_all.dta")


####
## Clean up the data
####

council_all<-subset(council_all, duplicate==F)
dim(council_all) # 8745

council_all$district[council_all$district=="atlarge" & council_all$state=="Texas"]<-"judge"
council_all$district[council_all$district=="Judge" & council_all$state=="Texas"]<-"judge"

council_all$district[council_all$district=="atlarge" & council_all$fips=="13135"]<-"Chairman"

council_all$district[council_all$district=="atlarge" & council_all$fips=="13151"]<- "Chairman"
council_all$district[council_all$district=="chair" & council_all$fips=="13151"]<- "Chairman"

council_all$district[council_all$district=="atlarge" & council_all$fips=="29099"]<- "presiding"

council_all<-subset(council_all, year>1985)

council_all <- group_by(council_all, fips, district) %>%
arrange(fips, year, district)%>%
    mutate(demshare.lag1 = lag(demshare, 1, along_with = year), 
     winner_pid.lag1 = lag(winner_pid, 1, along_with = year),
     winner_pid.lag2 = lag(winner_pid, 2, along_with = year),
	demshare.lag2 = lag(demshare, 2, along_with = year), 
            demshare.yearlag =lag(year, 1, along_with = year),
            demshare.yearlagD=as.numeric(year)-as.numeric(demshare.yearlag),
          demshare.Dlag1 = demshare.lag1-demshare
          )
council_all$demshare.lag1[council_all$ demshare.yearlagD>4]<-NA   


## Two Will County legislators per district after 2010, so elected every other cycle
## http://www.willcountyboard.com/about-the-board.html
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==17197 & council_all$year>2010]<-council_all$demshare.lag2[council_all$fips==17197& council_all$year>2010]
council_all$winner_pid.lag1[council_all$fips==17197& council_all$year>2010]<-council_all$winner_pid.lag2[council_all$fips==17197& council_all$year>2010]


## Two Champaign County legislators per district, so elected every other cycle
## http://www.co.champaign.il.us/countyboard
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==17019]<-council_all$demshare.lag2[council_all$fips==17019]
council_all$winner_pid.lag1[council_all$fips==17019]<-council_all$winner_pid.lag2[council_all$fips==17019]

## Two Dupage County legislators per district, so elected every other cycle
## https://www.dupageco.org/cobrd/
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==17043]<-council_all$demshare.lag2[council_all$fips==17043]
council_all$winner_pid.lag1[council_all$fips==17043]<-council_all$winner_pid.lag2[council_all$fips==17043]

## Orange County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37135 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37135 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37135 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37135 & council_all$district=="atlarge"]

## Union County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37179 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37179 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37179 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37179 & council_all$district=="atlarge"]

## Rowan County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37159 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37159 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37159 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37159 & council_all$district=="atlarge"]

## New Hanover County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37129 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37129 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37129 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37129 & council_all$district=="atlarge"]

## Iredell County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37097 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37097 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37097 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37097 & council_all$district=="atlarge"]

## Catawba County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37035 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37035 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37035 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37035 & council_all$district=="atlarge"]

## Cabarrus County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37025 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37025 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37025 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37025 & council_all$district=="atlarge"]

## Alamance County, NC at large legislators elected every other cycle
## So change lag1 to lag2 for demshare placebo check 
council_all$demshare.lag1[council_all$fips==37001 & council_all$district=="atlarge"]<-council_all$demshare.lag2[council_all$fips==37001 & council_all$district=="atlarge"]
council_all$winner_pid.lag1[council_all$fips==37001 & council_all$district=="atlarge"]<-council_all$winner_pid.lag2[council_all$fips==37001 & council_all$district=="atlarge"]

council_all$demseat<-NA
council_all$demseat[council_all$demshare.lag1>0]<- 1
council_all$demseat[council_all$demshare.lag1<0]<- 0

council_all$demseat2<-0
council_all$demseat2[council_all$demshare.lag1>0]<- 1
council_all$demseat2[council_all$demshare.lag1<0]<- -1

council_all$demseat3<-NA
council_all$demseat3[council_all$winner_pid.lag1== -1]<- 1
council_all$demseat3[council_all$winner_pid.lag1== 1]<- 0

pop<-read.dta( file="county_population_census.dta")

council_all2<-council_all
council_all2<-merge(council_all2, pop, by=c("fips", "year"), all.x=T)
pop<-subset(pop, year==2010)
pop<-dplyr::rename(pop, population_2010=population)  
pop<-dplyr::select(pop, -year) 
council_all2<-merge(council_all2, pop, by=c("fips"), all.x=T)

council_minyear <- council_all2 %>% group_by(fips) %>% slice(which.min(year)) %>% slice(which.max(population))
council_minyear<-rename(council_minyear, min_year_councilcands=year)
council_minyear<-subset(council_minyear,population>125000 )


# read in county fiscal data:
load(file="counties.Rdata")

data.panel<-data.use

data.panel <- group_by(data.panel, fips) %>%
    mutate(
    Total.Expenditure.PC.D2.log.lagged = lag.new(Total.Expenditure.PC.D2.log, 1, along_with = YearData),
    Total.Expenditure.PC.D3.log.lagged = lag.new(Total.Expenditure.PC.D3.log, 1, along_with = YearData),
     Total.Expenditure.PC.D2.lagged = lag.new(Total.Expenditure.PC.D2, 1, along_with = YearData),
          Total.Expenditure.PC.D3.lagged = lag.new(Total.Expenditure.PC.D3, 1, along_with = YearData),
     Total.Expenditure.PC.D1.lagged = lag.new(Total.Expenditure.PC.D1, 1, along_with = YearData))


## Council composition data. Has total number of seats in each council-year.
council_panel<-read.csv("councilcomp_export.csv") 
council_panel<-dplyr::select(council_panel, -population_est)
## this file has the number of seats in counties that are missing from our panel
council_seats<-read.csv("counties_missing_seats_total.csv")

council_panel<-merge(council_panel, council_seats, by.x=c( "fips"), by.y=c("fips"), all.x=T)

council_panel$seats_tot[is.na(council_panel$seats_tot)]<-council_panel$seats_tot2[is.na(council_panel$seats_tot)]

council_panel$partisan_diff<-abs(council_panel$seats_dem-council_panel$seats_rep)
##
## approximations
council_panel$partisan_diff[council_panel$fips==26163]<-13
council_panel$partisan_diff[council_panel$fips==17089]<-4

council_panel<-arrange(council_panel,fips, year)

council_panel <- group_by(council_panel, fips) %>%
    mutate(
            demshare_council2.lead1 = lead.new(demshare_council2, 1,along_with=year)     )    


years<-1990:2014
years<-as.data.frame(years)
colnames(years)<-"year"

fips<-unique(council_all$fips)
fips<-as.data.frame(fips)
colnames(fips)<-"fips"
fips<-merge(fips, years, all.x=T, all.y=T)
council_panel<-merge(council_panel,fips, by=c("fips", "year"), all.x=T, all.y=T)

council_panel <- group_by(council_panel, fips) %>%
    mutate(
    partisan_diff.L1 = lead.new(partisan_diff, 1, along_with = year),
      seats_dem.L1 = lead.new(seats_dem, 1, along_with = year),
        seats_rep.L1 = lead.new(seats_rep, 1, along_with = year),
          seats_tot.L1 = lead.new(seats_tot, 1, along_with = year))

   

data.panel<-merge(data.panel, council_panel, by.x=c("YearData", "fips"), by.y=c("year", "fips"))
data.panel<-merge(data.panel, pop, by.x=c("fips"), by.y=c( "fips"))

data.panel <- group_by(data.panel, fips) %>%
    mutate(
    demshare_council.lag1 = lag.new(demshare_council2, 1, along_with = YearData), 
    demshare_council.D=demshare_council2-demshare_council.lag1)


data.panel$partisan_diff2<-(data.panel$seats_dem-data.panel$seats_rep)



data.panel$state_year<-factor(paste(data.panel$STATEAB, data.panel$YearData))

data.panel$year<-factor(data.panel$YearData)

data.panel<-subset(data.panel, YearData>1989)
data.panel<-subset(data.panel, population_2010>150000)

data.panel<-subset(data.panel, STATEAB!="CA" )
data.panel<-subset(data.panel, STATEAB!="MN" )
data.panel<-subset(data.panel, STATEAB!="WI" )
data.panel<-subset(data.panel, STATEAB!="AK" )
data.panel<-subset(data.panel, STATEAB!="SD" )
data.panel<-subset(data.panel, STATEAB!="ND" )
data.panel<-subset(data.panel, STATEAB!="LA" )
data.panel<-subset(data.panel, fips!="12073" )
data.panel<-subset(data.panel, fips!="12127" )
## lane county oregon is nonpartisan
data.panel<-subset(data.panel, fips!="41039" )

####
## Figure 6: Association between Changes in Partisan Composition of County Legislators and Changes in County Spending. 
#####
width<-.0333
data.panel <- mutate(data.panel, bin=cut(demshare_council.D, breaks=seq(-1,1, width)))
bin1.df <- data.frame(bin.mean=tapply(data.panel$Total.Expenditure.PC.D2.log.lagged,
                         data.panel$bin, mean, na.rm=TRUE),
                     mid=seq(-1 + width/2, 1 - width/2, width),
                     n=tapply(data.panel$Total.Expenditure.PC.D2.log.lagged, data.panel$bin, length))
                     
c <- ggplot(data.panel, aes(demshare_council.D, Total.Expenditure.PC.D2.log.lagged))
c <- c + geom_point(data=bin1.df, aes(x=mid, y=bin.mean, size=log(n)), pch=1)
c <- c + stat_smooth(method = "lm", formula = y ~ x , size = 1, color="black")
 c <- c + theme_bw() 
c <- c + scale_x_continuous(breaks=c(-.75,-.5,-.25,0,.25,.5,.75), limits=c(-.75,.75))
c <- c +  coord_cartesian(ylim=c(-.08,.08),xlim=c(-.8,.7))
c <- c +labs(size = 'Logged # of Obs.\nin 3.3% Bin',
        x = "Change in Democratic Seat Share",
        y = "Change in Logged Spending Per Capita")
pdf(file="figures/figure6_policy_firstdifferences_logged.pdf", height=4, width=7.5)
c
dev.off()



data.use <- group_by(data.use, govid) %>%
    mutate(          
    	Total.Expenditure.PC.lagged.D1 =Total.Expenditure.PC- Total.Expenditure.PC.Lag1,
    	Total.Expenditure.PC.lagged.D1.log =log(Total.Expenditure.PC+.01)- log(Total.Expenditure.PC.Lag1+.01))
           
states<-read.csv(file="states_icpsr.csv",
                 stringsAsFactors=FALSE)

states<-states[,c(1,4)]

icma_institutions_imputed<-read.csv(file="countyfog_imputedfully_long.csv",
                                    stringsAsFactors=FALSE)
## Do some manipulations to make full fips codes for ICMA data
icma_institutions_imputed<-rename(icma_institutions_imputed, abb=STATE)
icma_institutions_imputed<-merge(icma_institutions_imputed, states,
                                 by.x="abb", by.y="STATEAB")
                                 
icma_institutions_imputed$fips<-NA
icma_institutions_imputed$fips[
        nchar(icma_institutions_imputed$FIPS_COUNTY) ==2 ]<-
             paste( icma_institutions_imputed$fips_state[nchar(icma_institutions_imputed$FIPS_COUNTY)==2], "0",
                   icma_institutions_imputed$FIPS_COUNTY[nchar(icma_institutions_imputed$FIPS_COUNTY)==2], sep="")

icma_institutions_imputed$fips[nchar(icma_institutions_imputed$FIPS_COUNTY)==3] <-
    paste(icma_institutions_imputed$fips_state[nchar(icma_institutions_imputed$FIPS_COUNTY)==3],
          icma_institutions_imputed$FIPS_COUNTY[nchar(icma_institutions_imputed$FIPS_COUNTY)==3], sep="")
          
icma_institutions_imputed$fips[nchar(icma_institutions_imputed$FIPS_COUNTY)==1] <-
    paste(icma_institutions_imputed$fips_state[nchar(icma_institutions_imputed$FIPS_COUNTY)==1], "00",
          icma_institutions_imputed$FIPS_COUNTY[nchar(icma_institutions_imputed$FIPS_COUNTY)==1], sep="")
icma_institutions_imputed$fips<-as.numeric(as.vector(icma_institutions_imputed$fips))

colnames(icma_institutions_imputed)[colnames(icma_institutions_imputed)=="fog"]<-"icma_fog"
colnames(icma_institutions_imputed)[colnames(icma_institutions_imputed)=="initiative"]<-"icma_initiative"
colnames(icma_institutions_imputed)[colnames(icma_institutions_imputed)=="partisan"]<-"icma_partisan"
colnames(icma_institutions_imputed)[colnames(icma_institutions_imputed)=="referendum"]<-"icma_referendum"
icma_institutions_imputed<-dplyr::select(icma_institutions_imputed, -abb)



# merge fiscal data with council elections data
dim(council_all)
council_all$id<-1:dim(council_all)[1]

data_rd<-merge(data.use, council_all, by.x=c("YearData", "fips"), by.y=c("year", "fips"))
dim(data_rd)
data_rd<-merge(data_rd, council_panel, by.x=c("YearData", "fips"), by.y=c("year", "fips"), all.x=T)
dim(data_rd)

council_seats<-read.csv("counties_missing_seats_total.csv")
council_seats<-rename(council_seats, seats_tot3=seats_tot2)
data_rd<-merge(data_rd, council_seats, by.x=c( "fips"), by.y=c( "fips"), all.x=T)
data_rd$seats_tot[is.na(data_rd$seats_tot)]<-data_rd$seats_tot3[is.na(data_rd$seats_tot)]

icma_institutions_imputed$fips<-as.numeric(icma_institutions_imputed$fips)

data_rd<-merge(data_rd, icma_institutions_imputed, by.x=c("YearData", "fips"), by.y=c("year", "fips"), all.x=T)
dim(data_rd)

icma_fog_corrections<-read.csv(file="icma_fog_corrections.csv")


## import and merge the 1987 CoG for FoG data
counties_1987census<-read.dta(file="counties_1987census.dta")
data_rd<-merge(data_rd, counties_1987census, by="fips", all.x=T)
data_rd<-merge(data_rd, icma_fog_corrections, by="fips", all.x=T)

data_rd$fog_combined<-data_rd$icma_fog
data_rd$fog_combined[is.na(data_rd$fog_combined)]<-data_rd$fog[is.na(data_rd$fog_combined)]
data_rd$fog_combined[!is.na(data_rd$fog_update)]<-data_rd$fog_update[!is.na(data_rd$fog_update)]

data_rd$fog_combined[data_rd$fips=="17161"]<-2

### Import and merge information on county fiscal responsibilities
#states_county_expenditure_areas<-read.csv(file="states_county_expenditure_areas.csv")
#data_rd<-merge(data_rd, states_county_expenditure_areas, by.x="STATEAB", by.y="abb", all.x=T)
#dim(data_rd)


#data_rd$duplicate<-duplicated(paste(data_rd$govid, data_rd$YearData, data_rd$demshare))

data_rd_popest<-subset(data_rd, YearData==2010)
data_rd_popest<-group_by(data_rd_popest, fips)
data_rd_popest<-dplyr::summarize(data_rd_popest, population_2010_est=mean(population_est))  

data_rd_popest<-dplyr::select(data_rd_popest, fips, population_2010_est)
data_rd<-merge(data_rd, pop, by="fips", all.x=T)
data_rd<-merge(data_rd, data_rd_popest, by="fips", all.x=T)
data_rd$population_2010[is.na(data_rd$population_2010)]<-data_rd$population_2010_est[is.na(data_rd$population_2010)]


council_size <- data_rd %>% group_by(fips) %>% slice(which.max(seats_tot))
council_size<-dplyr::select(council_size, fips, seats_tot)
council_size<-dplyr::rename(council_size, seats_tot_max=seats_tot)   
data_rd<-merge(data_rd, council_size, by="fips", all.x=T)
data_rd$seats_tot[is.na(data_rd$seats_tot)]<-data_rd$seats_tot_max[is.na(data_rd$seats_tot)]

## Subset to data we use for analysis
#data_rd<-subset(data_rd, population_est>150000 )
dim(data_rd)
data_rd<-subset(data_rd, YearData>1989 )
data_rd<-subset(data_rd, YearData<2014 )
dim(data_rd)
#data_rd<-subset(data_rd,!is.na(demshare))
#dim(data_rd)
data_rd$cluster<-paste(data_rd$fips, data_rd$YearData)


###
## drop states with all nonpartisan elections and clean up the data
###
data_rd<-subset(data_rd, STATEAB!="CA" )
data_rd<-subset(data_rd, STATEAB!="MN" )
data_rd<-subset(data_rd, STATEAB!="WI" )
data_rd<-subset(data_rd, STATEAB!="AK" )
data_rd<-subset(data_rd, STATEAB!="SD" )
data_rd<-subset(data_rd, STATEAB!="ND" )
data_rd<-subset(data_rd, STATEAB!="LA" )
data_rd<-subset(data_rd, fips!="12073" )
data_rd<-subset(data_rd, fips!="12127" )
## lane county oregon is nonpartisan
data_rd<-subset(data_rd, fips!="41039" )
## washington county oregon is nonpartisan
data_rd<-subset(data_rd, fips!="41067" )
### NEW
## lane county oregon is nonpartisan
data_rd<-subset(data_rd, fips!="41039" )
## Bibb county GA is consolidated
data_rd<-subset(data_rd, fips!="13021" )
## Whatcom County, Washington is nonpartisan
data_rd<-subset(data_rd, fips!="53073" )

## Make sure there are no duplicate rows in the data
data_rd$duplicate<-duplicated(paste(data_rd$fips, data_rd$YearData, data_rd$district))
data_rd<-subset(data_rd, duplicate==F)
dim(data_rd) # 6070


## Drop some screwy counties
## Calvert County MD has very strange rules before 2003
drop<-data_rd$fips==24009 & data_rd$YearData<2003
data_rd<-subset(data_rd, drop==F)

## create variables that average the change in fiscal policies in t+2 and t+3 after elections
data_rd<-group_by(data_rd, fips) %>%
	mutate( 
	Total.Expenditure.PC.D23.log.old=(Total.Expenditure.PC.D2.log+Total.Expenditure.PC.D3.log)/2,
	Total.Expenditure.PC.D23.log=rowMeans(cbind(Total.Expenditure.PC.D2.log,Total.Expenditure.PC.D3.log), na.rm=T),
	Total.Allocational.PC.D23.log=rowMeans(cbind(Total.Allocational.PC.D2.log,Total.Allocational.PC.D3.log), na.rm=T),
	Total.PublicProtection.PC.D23.log=rowMeans(cbind(Total.PublicProtection.PC.D2.log,Total.PublicProtection.PC.D3.log), na.rm=T),
	Total.AdminMisc.PC.D23.log=rowMeans(cbind(Total.AdminMisc.PC.D2.log,Total.AdminMisc.PC.D3.log), na.rm=T),
	Total.ParksResources.PC.D23.log=rowMeans(cbind(Total.ParksResources.PC.D2.log,Total.ParksResources.PC.D3.log), na.rm=T),
	Total.RoadsParking.PC.D23.log=rowMeans(cbind(Total.RoadsParking.PC.D2.log,Total.RoadsParking.PC.D3.log), na.rm=T),
	Total.SanitationUtilities.PC.D23.log=rowMeans(cbind(Total.SanitationUtilities.PC.D2.log,Total.SanitationUtilities.PC.D3.log), na.rm=T),
	Total.EducationLibraries.PC.D23.log=rowMeans(cbind(Total.EducationLibraries.PC.D2.log,Total.EducationLibraries.PC.D3.log), na.rm=T),
	Total.Developmental.PC.D23.log=rowMeans(cbind(Total.Developmental.PC.D2.log,Total.Developmental.PC.D3.log), na.rm=T),
	Total.Redistribution.PC.D23.log=rowMeans(cbind(Total.Redistribution.PC.D2.log,Total.Redistribution.PC.D3.log), na.rm=T),
	Total.Expenditure.real.D23.log=rowMeans(cbind(Total.Expenditure.real.D2.log,Total.Expenditure.real.D3.log), na.rm=T),
	Total.Expenditure.PC.D23=rowMeans(cbind(Total.Expenditure.PC.D2,Total.Expenditure.PC.D3), na.rm=T),
	Total.Expenditure.PC.D234.log=rowMeans(cbind(Total.Expenditure.PC.D2.log,Total.Expenditure.PC.D3.log,Total.Expenditure.PC.D4.log), na.rm=T),
	Total.Expenditure.PC.D23.log.col=rowMeans(cbind(Total.Expenditure.PC.D2.log.col,Total.Expenditure.PC.D3.log.col), na.rm=T),
	corrections.PC.D23.log=rowMeans(cbind(corrections.PC.D2.log,corrections.PC.D3.log), na.rm=T),
	naturalresources.PC.D23.log=rowMeans(cbind(naturalresources.PC.D2.log,naturalresources.PC.D3.log), na.rm=T),
	nec.PC.D23.log=rowMeans(cbind(nec.PC.D2.log,nec.PC.D3.log), na.rm=T),
	charges.PC.D23.log=rowMeans(cbind(charges.PC.D2.log,charges.PC.D3.log), na.rm=T),
	revenueig.PC.D23.log=rowMeans(cbind(revenueig.PC.D2.log,revenueig.PC.D3.log), na.rm=T),
	education.PC.D23.log=rowMeans(cbind(education.PC.D2.log,education.PC.D3.log), na.rm=T),
	fire.PC.D23.log=rowMeans(cbind(fire.PC.D2.log,fire.PC.D3.log), na.rm=T),
	police.PC.D23.log=rowMeans(cbind(police.PC.D2.log,police.PC.D3.log), na.rm=T),
	health.PC.D23.log=rowMeans(cbind(health.PC.D2.log,health.PC.D3.log), na.rm=T),
	highways.PC.D23.log=rowMeans(cbind(highways.PC.D2.log,highways.PC.D3.log), na.rm=T),
	housing.PC.D23.log=rowMeans(cbind(housing.PC.D2.log,housing.PC.D3.log), na.rm=T),
	libraries.PC.D23.log=rowMeans(cbind(libraries.PC.D2.log,libraries.PC.D3.log), na.rm=T),
	parks.PC.D23.log=rowMeans(cbind(parks.PC.D2.log,parks.PC.D3.log), na.rm=T),
	sanitation.PC.D23.log=rowMeans(cbind(sanitation.PC.D2.log,sanitation.PC.D3.log), na.rm=T),
	utilities.PC.D23.log=rowMeans(cbind(utilities.PC.D2.log,utilities.PC.D3.log), na.rm=T),
	welfare.PC.D23.log=rowMeans(cbind(welfare.PC.D2.log,welfare.PC.D3.log), na.rm=T),
	interest.PC.D23.log=rowMeans(cbind(interest.PC.D2.log,interest.PC.D3.log), na.rm=T),
	admin.PC.D23.log=rowMeans(cbind(admin.PC.D2.log,admin.PC.D3.log), na.rm=T),
	revenue.PC.D23.log=rowMeans(cbind(revenue.PC.D2.log,revenue.PC.D3.log), na.rm=T),
	revenueown.PC.D23.log=rowMeans(cbind(revenueown.PC.D2.log,revenueown.PC.D3.log), na.rm=T),
	Total.Taxes.D23.log=rowMeans(cbind(Total.Taxes.D2.log,Total.Taxes.D3.log), na.rm=T),
	salestax.PC.D23.log=rowMeans(cbind(salestax.PC.D2.log,salestax.PC.D3.log), na.rm=T),
	salestax.select.PC.D23.log=rowMeans(cbind(salestax.select.PC.D2.log,salestax.select.PC.D3.log), na.rm=T),
	salestax.general.PC.D23.log=rowMeans(cbind(salestax.general.PC.D2.log,salestax.general.PC.D3.log), na.rm=T),
	propertytax.PC.D23.log=rowMeans(cbind(propertytax.PC.D2.log,propertytax.PC.D3.log), na.rm=T),
	debt.PC.D23.log=rowMeans(cbind(debt.PC.D2.log,debt.PC.D3.log), na.rm=T)
)

data_rd_100kthreshold<-subset(data_rd, population_2010>100000 )
data_rd<-subset(data_rd, population_2010>150000 )
## chart of size of county legislators
data_councilsize<-group_by(data_rd, fips) %>%
	summarise(seats=mean(seats_tot))
mean_countycouncil_size<-mean(data_councilsize$seats, na.rm=T)

	
####
## create weights for the size of the county council
####
data_rd$mean_seats_tot<-mean(data_rd$seats_tot, na.rm=T)
data_rd$weight<-mean_countycouncil_size/data_rd$seats_tot
data_rd<-group_by(data_rd, fips, YearData) %>%
	mutate(sum_weights=sum(weight))

data_rd_100kthreshold$mean_seats_tot<-mean(data_rd_100kthreshold$seats_tot, na.rm=T)
data_rd_100kthreshold$weight<-mean_countycouncil_size/data_rd_100kthreshold$seats_tot
data_rd_100kthreshold<-group_by(data_rd_100kthreshold, fips, YearData) %>%
	mutate(sum_weights=sum(weight))

dim(data_rd)
length(unique(data_rd$fips))

####
##Figure 1: Size of County Legislatures in our Dataset
####
pdf("figures/figure1_Council_Size.pdf", width=7.5, height=4)
(ggplot(data_councilsize, aes(seats)) +
		geom_histogram(binwidth = 1)+
		scale_y_continuous("Number of Counties")+
		scale_x_continuous("Size of Legislature",breaks=seq(0,40,5))+
		theme_bw())
dev.off()
        
#pdf("figures/Data_LegSize.pdf", width=8, height=4.5)
#(ggplot(data_rd, aes(seats_tot)) +
# geom_histogram()+
#   scale_y_continuous("Number of Observations in Dataset")+
#    scale_x_continuous("Size of Legislature",breaks=seq(0,40,5))+
#        theme_bw())
 #       dev.off()


# export for descriptives:
county_desc_summ <- ddply(data_rd, ~ fips, summarize,
      yearmax = max(YearData),
      yearmin = min(YearData))
write.csv(county_desc_summ,"counties_map.csv",row.names=F)


source("rd.export.2017.R")


####
## Effect of council partisanship on expenditures
## There are several tables and figures outputted. See below.
####

cct.council.exp.log <- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))


summary(with(data_rd,Total.Expenditure.PC))

summary(with(data_rd,education.PC+fire.PC+police.PC+health.PC+highways.PC+housing.PC+libraries.PC+parks.PC+sanitation.PC+utilities.PC+welfare.PC+interest.PC+admin.PC+corrections.PC+naturalresources.PC+nec.PC))


summary(with(data_rd, rdrobust(y=(education.PC.D23.log+fire.PC.D23.log+police.PC.D23.log+health.PC.D23.log+highways.PC.D23.log+housing.PC.D23.log+libraries.PC.D23.log+parks.PC.D23.log+sanitation.PC.D23.log+utilities.PC.D23.log+welfare.PC.D23.log+interest.PC.D23.log+admin.PC.D23.log+corrections.PC.D23.log+naturalresources.PC.D23.log+nec.PC.D23.log)/16 , x=demshare, cluster=cluster, weight=weight)))

cct.education.log <- with(data_rd, rdrobust(y=education.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.nec.log <- with(data_rd, rdrobust(y=nec.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.corrections.log <- with(data_rd, rdrobust(y=corrections.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.naturalresources.log <- with(data_rd, rdrobust(y=naturalresources.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))


cct.fire.log <- with(data_rd, rdrobust(y=fire.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.police.log <- with(data_rd, rdrobust(y=police.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.health.log <- with(data_rd, rdrobust(y=health.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.highways.log <- with(data_rd, rdrobust(y=highways.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.housing.log <- with(data_rd, rdrobust(y=housing.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.libraries.log <- with(data_rd, rdrobust(y=libraries.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.parks.log <- with(data_rd, rdrobust(y=parks.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.sanitation.log <- with(data_rd, rdrobust(y=sanitation.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.utilities.log <- with(data_rd, rdrobust(y=utilities.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.welfare.log <- with(data_rd, rdrobust(y=welfare.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.interest.log <- with(data_rd, rdrobust(y=interest.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.admin.log <- with(data_rd, rdrobust(y=admin.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.council.exp.log2<- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log  , x=demshare, cluster=cluster, weight=weight, level=90))

cct.education.log2 <- with(data_rd, rdrobust(y=education.PC.D23.log  , x=demshare, cluster=cluster, weight=weight , level=90))

cct.nec.log2 <- with(data_rd, rdrobust(y=nec.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.corrections.log2 <- with(data_rd, rdrobust(y=corrections.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.naturalresources.log2 <- with(data_rd, rdrobust(y=naturalresources.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.fire.log2 <- with(data_rd, rdrobust(y=fire.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.police.log2 <- with(data_rd, rdrobust(y=police.PC.D23.log, x=demshare, cluster=cluster, weight=weight , level=90))

cct.health.log2 <- with(data_rd, rdrobust(y=health.PC.D23.log , x=demshare, cluster=cluster, weight=weight , level=90))

cct.highways.log2 <- with(data_rd, rdrobust(y=highways.PC.D23.log, x=demshare, cluster=cluster, weight=weight , level=90))

cct.housing.log2 <- with(data_rd, rdrobust(y=housing.PC.D23.log , x=demshare, cluster=cluster, weight=weight , level=90))

cct.libraries.log2 <- with(data_rd, rdrobust(y=libraries.PC.D23.log , x=demshare, cluster=cluster, weight=weight , level=90))

cct.parks.log2 <- with(data_rd, rdrobust(y=parks.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.sanitation.log2 <- with(data_rd, rdrobust(y=sanitation.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.utilities.log2 <- with(data_rd, rdrobust(y=utilities.PC.D23.log , x=demshare, cluster=cluster, weight=weight , level=90))

cct.welfare.log2 <- with(data_rd, rdrobust(y=welfare.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.interest.log2 <- with(data_rd, rdrobust(y=interest.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.admin.log2 <- with(data_rd, rdrobust(y=admin.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90))

coefs = data.frame(rbind(cbind(cct.council.exp.log$coef[1,1],cct.council.exp.log$ci[3,1],cct.council.exp.log$ci[3,2],cct.council.exp.log2$ci[3,1],cct.council.exp.log2$ci[3,2]),
                         #cbind(cct.expenditures.general.log$coef[1,1],cct.expenditures.general.log$ci[3,1],cct.expenditures.general.log$ci[3,2]),
               
                         cbind(cct.education.log$coef[1,1],cct.education.log$ci[3,1],cct.education.log$ci[3,2],cct.education.log2$ci[3,1],cct.education.log2$ci[3,2]),
                         cbind(cct.fire.log$coef[1,1],cct.fire.log$ci[3,1],cct.fire.log$ci[3,2],cct.fire.log2$ci[3,1],cct.fire.log2$ci[3,2]),
                         cbind(cct.police.log$coef[1,1],cct.police.log$ci[3,1],cct.police.log$ci[3,2],cct.police.log2$ci[3,1],cct.police.log2$ci[3,2]),
                                        cbind(cct.corrections.log$coef[1,1],cct.corrections.log$ci[3,1],cct.corrections.log$ci[3,2],cct.corrections.log2$ci[3,1],cct.corrections.log2$ci[3,2]),
                 cbind(cct.health.log$coef[1,1],cct.health.log$ci[3,1],cct.health.log$ci[3,2],cct.health.log2$ci[3,1],cct.health.log2$ci[3,2]),
                         cbind(cct.highways.log$coef[1,1],cct.highways.log$ci[3,1],cct.highways.log$ci[3,2],cct.highways.log2$ci[3,1],cct.highways.log2$ci[3,2]),
                         cbind(cct.housing.log$coef[1,1],cct.housing.log$ci[3,1],cct.housing.log$ci[3,2],cct.housing.log2$ci[3,1],cct.housing.log2$ci[3,2]),
                         cbind(cct.libraries.log$coef[1,1],cct.libraries.log$ci[3,1],cct.libraries.log$ci[3,2],cct.libraries.log2$ci[3,1],cct.libraries.log2$ci[3,2]),
                         cbind(cct.parks.log$coef[1,1],cct.parks.log$ci[3,1],cct.parks.log$ci[3,2],cct.parks.log2$ci[3,1],cct.parks.log2$ci[3,2]),
                                         cbind(cct.naturalresources.log$coef[1,1],cct.naturalresources.log$ci[3,1],cct.naturalresources.log$ci[3,2],cct.naturalresources.log2$ci[3,1],cct.naturalresources.log2$ci[3,2]),
        cbind(cct.sanitation.log$coef[1,1],cct.sanitation.log$ci[3,1],cct.sanitation.log$ci[3,2],cct.sanitation.log2$ci[3,1],cct.sanitation.log2$ci[3,2]),
                         cbind(cct.utilities.log$coef[1,1],cct.utilities.log$ci[3,1],cct.utilities.log$ci[3,2],cct.utilities.log2$ci[3,1],cct.utilities.log2$ci[3,2]),
                         cbind(cct.welfare.log$coef[1,1],cct.welfare.log$ci[3,1],cct.welfare.log$ci[3,2],cct.welfare.log2$ci[3,1],cct.welfare.log2$ci[3,2]),
                         cbind(cct.interest.log$coef[1,1],cct.interest.log$ci[3,1],cct.interest.log$ci[3,2],cct.interest.log2$ci[3,1],cct.interest.log2$ci[3,2]),
                         cbind(cct.admin.log$coef[1,1],cct.admin.log$ci[3,1],cct.admin.log$ci[3,2],cct.admin.log2$ci[3,1],cct.admin.log2$ci[3,2]),
                         cbind(cct.nec.log$coef[1,1],cct.nec.log$ci[3,1],cct.nec.log$ci[3,2],cct.nec.log2$ci[3,1],cct.nec.log2$ci[3,2])
))
names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("Total Expenditures",
     #          "General Expenditures",
               "Education",
               "Fire",
               "Police",
                "Corrections",
              "Health",
               "Highways",
               "Housing",
               "Libraries",
               "Parks",
                "Natural Resources",
              "Sanitation",
               "Utilities",
               "Welfare",
               "Interest",
               "Administration",
               "Misc.")
coefs[2:17,]<-coefs[2:17,][order(-coefs$coef[2:17]),]

####
## Effect of council partisanship on all expenditure categories (extra figure)
####

pdf("figures/extrafigure_coefplot_exp_log_allcategories.pdf",height=5,width=7)
ggplot(coefs, aes(c(17:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
      geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                     colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient + Robust CI\n(Ave. in Years 2 and 3 after Election)",breaks=seq(-0.4, .8, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  max(coefs$ci95.hi)),
                         labels=seq(-0.4, .8, 0.1)) +
      scale_x_continuous("Outcome Variable (logged per capita)",breaks=c(17:1),
                         limits=c(1,17.5),
                         labels=coefs$vars) + 
      coord_flip() + 
      geom_text(nudge_x=.5, size=2) +
      theme_bw()
dev.off()



cct.Allocational.log <- with(data_rd, rdrobust(y= Total.Allocational.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
cct.Developmental.log <- with(data_rd, rdrobust(y= Total.Developmental.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
cct.Redistribution.log <- with(data_rd, rdrobust(y= Total.Redistribution.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

cct.Allocational.log2 <- with(data_rd, rdrobust(y= Total.Allocational.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90))
cct.Developmental.log2 <- with(data_rd, rdrobust(y= Total.Developmental.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90))
cct.Redistribution.log2 <- with(data_rd, rdrobust(y= Total.Redistribution.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90))


coefs = data.frame(rbind(cbind(cct.council.exp.log$coef[1,1],cct.council.exp.log$ci[3,1],cct.council.exp.log$ci[3,2],cct.council.exp.log2$ci[3,1],cct.council.exp.log2$ci[3,2]),
cbind(cct.Allocational.log$coef[1,1],cct.Allocational.log$ci[3,1],cct.Allocational.log$ci[3,2],cct.Allocational.log2$ci[3,1],cct.Allocational.log2$ci[3,2]),
                         #cbind(cct.expenditures.general.log$coef[1,1],cct.expenditures.general.log$ci[3,1],cct.expenditures.general.log$ci[3,2]),
               
                         cbind(cct.Developmental.log$coef[1,1],cct.Developmental.log$ci[3,1],cct.Developmental.log$ci[3,2],cct.Developmental.log2$ci[3,1],cct.Developmental.log2$ci[3,2]),
   cbind(cct.education.log$coef[1,1],cct.education.log$ci[3,1],cct.education.log$ci[3,2],cct.education.log2$ci[3,1],cct.education.log2$ci[3,2]),
                         cbind(cct.Redistribution.log$coef[1,1],cct.Redistribution.log$ci[3,1],cct.Redistribution.log$ci[3,2],cct.Redistribution.log2$ci[3,1],cct.Redistribution.log2$ci[3,2])))

names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("Total Expenditures",
     #          "General Expenditures",
               "Allocational
(fire, police, 
corrections, interest, 
administrative, parks, 
natural resources, 
libraries, sanitation, misc)",
               "Developmental
               (highways, parking, 
utilities)",
               "Education",
               "Redistribution
(health, hospitals, welfare, housing")

####
## Figure D1: Aggregated spending categories from Peterson (1981)
####

pdf("figures/figureD1_coefplot_exp_log_petersoncategories.pdf",height=5,width=7)
ggplot(coefs, aes(c(5:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
      geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                     colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient + Robust CI\n(Ave. in Years 2 and 3 after Election)",breaks=seq(-0.4, .8, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  max(coefs$ci95.hi)),
                         labels=seq(-0.4, .8, 0.1)) +
      scale_x_continuous("Outcome Variable (logged per capita)",breaks=c(5:1),
                         limits=c(1,5.2),
                         labels=coefs$vars,position = "bottom") + 
      coord_flip() + 
      geom_text(nudge_x=.15, size=4) +
      theme_bw()
dev.off()

## Expenditures, logged, D2+D3
cct.council.exp.log <- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))

cct.PublicProtection.log <- with(data_rd, rdrobust(y=Total.PublicProtection.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.Redistribution.log <- with(data_rd, rdrobust(y=Total.Redistribution.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.education.log <- with(data_rd, rdrobust(y=Total.EducationLibraries.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.highways.log <- with(data_rd, rdrobust(y=highways.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.ParksResources.log <- with(data_rd, rdrobust(y=Total.ParksResources.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.sanitation.log <- with(data_rd, rdrobust(y=Total.SanitationUtilities.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.interest.log <- with(data_rd, rdrobust(y=interest.PC.D23.log, x=demshare, cluster=cluster, weight=weight ))

cct.AdminMisc.log <- with(data_rd, rdrobust(y=Total.AdminMisc.PC.D23.log , x=demshare, cluster=cluster, weight=weight ))

cct.council.exp.log2 <- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90))

cct.PublicProtection.log2 <- with(data_rd, rdrobust(y=Total.PublicProtection.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.Redistribution.log2 <- with(data_rd, rdrobust(y=Total.Redistribution.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.education.log2 <- with(data_rd, rdrobust(y=Total.EducationLibraries.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.highways.log2 <- with(data_rd, rdrobust(y=highways.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.ParksResources.log2 <- with(data_rd, rdrobust(y=Total.ParksResources.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.sanitation.log2 <- with(data_rd, rdrobust(y=Total.SanitationUtilities.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.interest.log2 <- with(data_rd, rdrobust(y=interest.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90 ))

cct.AdminMisc.log2 <- with(data_rd, rdrobust(y=Total.AdminMisc.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90 ))



coefs = data.frame(rbind(cbind(cct.council.exp.log$coef[1,1],cct.council.exp.log$ci[3,1],cct.council.exp.log$ci[3,2],cct.council.exp.log2$ci[3,1],cct.council.exp.log2$ci[3,2]),

cbind(cct.PublicProtection.log$coef[1,1],cct.PublicProtection.log$ci[3,1],cct.PublicProtection.log$ci[3,2],cct.PublicProtection.log2$ci[3,1],cct.PublicProtection.log2$ci[3,2]),
cbind(cct.Redistribution.log$coef[1,1],cct.Redistribution.log$ci[3,1],cct.Redistribution.log$ci[3,2],cct.Redistribution.log2$ci[3,1],cct.Redistribution.log2$ci[3,2]),
cbind(cct.education.log$coef[1,1],cct.education.log$ci[3,1],cct.education.log$ci[3,2],cct.education.log2$ci[3,1],cct.education.log2$ci[3,2]),
cbind(cct.highways.log$coef[1,1],cct.highways.log$ci[3,1],cct.highways.log$ci[3,2],cct.highways.log2$ci[3,1],cct.highways.log2$ci[3,2]),
cbind(cct.ParksResources.log$coef[1,1],cct.ParksResources.log$ci[3,1],cct.ParksResources.log$ci[3,2],cct.ParksResources.log2$ci[3,1],cct.ParksResources.log2$ci[3,2]),
cbind(cct.sanitation.log$coef[1,1],cct.sanitation.log$ci[3,1],cct.sanitation.log$ci[3,2],cct.sanitation.log2$ci[3,1],cct.sanitation.log2$ci[3,2]),
cbind(cct.interest.log$coef[1,1],cct.interest.log$ci[3,1],cct.interest.log$ci[3,2],cct.interest.log2$ci[3,1],cct.interest.log2$ci[3,2]),
cbind(cct.AdminMisc.log$coef[1,1],cct.AdminMisc.log$ci[3,1],cct.AdminMisc.log$ci[3,2],cct.AdminMisc.log2$ci[3,1],cct.AdminMisc.log2$ci[3,2])))

names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("Total Expenditures",
				"Police, Fire, and Corrections",
				"Redistribution (Health, Hospitals, Housing, & Welfare)",
				"Education and Libraries",
				"Roads",
				"Parks and Natural Resources",
				"Sanitation and Utilities",
				"Interest",
				"Admin. and Misc.")
				
coefs[2:9,]<-coefs[2:9,][order(-coefs$coef[2:9]),]

####
## Table D1: Expenditure Results
####
rd.export(list(cct.council.exp.log,cct.PublicProtection.log,cct.Redistribution.log,cct.education.log,cct.highways.log,cct.ParksResources.log,cct.sanitation.log,cct.interest.log,cct.AdminMisc.log))


####
## Figure 9: The effect of legislative partisanship on changes in logged program-area per capita expenditures in the fiscal years two and three years after an election.
####

pdf("figures/figure9_coefplot_exp_log_collapsedcategories.pdf",height=5,width=7.5)
ggplot(coefs, aes(c(9:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
      geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                     colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient + Robust CI\n(Ave. in Years 2 and 3 after Election)",breaks=seq(-0.4, .8, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  max(coefs$ci95.hi)),
                         labels=seq(-0.4, .8, 0.1)) +
      scale_x_continuous("Outcome Variable (logged per capita)",breaks=c(9:1),
                         limits=c(1,9.3),
                         labels=coefs$vars,position = "bottom") + 
      coord_flip() + 
       geom_text(nudge_x=.3, size=3) +
     theme_bw()
dev.off()



####
## Effect of council partisanship on revenues
## There are several tables and figures outputted. See below.
####

cct.council.revenues.log <- with(data_rd, rdrobust(y=revenue.PC.D23.log , x=demshare, cluster=cluster, weight=weight))

summary(with(data_rd, rdrobust(y=(Total.Taxes.D23.log+charges.PC.D23.log+revenueig.PC.D23.log)/3 , x=demshare, cluster=cluster, weight=weight)))

cct.council.revenuesown.log <- with(data_rd, rdrobust(y=revenueown.PC.D23.log , x=demshare, cluster=cluster, weight=weight))

 cct.council.taxes.log <- with(data_rd, rdrobust(y=Total.Taxes.D23.log, x=demshare, cluster=cluster, weight=weight))
	
 cct.council.salestaxes.log <- with(data_rd, rdrobust(y=salestax.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
 
  cct.council.salestaxes.select.log <- with(data_rd, rdrobust(y=salestax.select.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
  
 cct.council.salestaxes.general.log <- with(data_rd, rdrobust(y=salestax.general.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

 cct.council.propertytaxes.log <- with(data_rd, rdrobust(y=propertytax.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
 
  cct.council.debt.log <- with(data_rd, rdrobust(y= debt.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

  cct.council.charges.log <- with(data_rd, rdrobust(y= charges.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
  cct.council.revenueig.log <- with(data_rd, rdrobust(y= revenueig.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
 
####
## Table D2: Revenue Results.
####

rd.export(list(cct.council.revenues.log,cct.council.revenuesown.log, cct.council.taxes.log,cct.council.salestaxes.log, cct.council.propertytaxes.log , cct.council.charges.log, cct.council.revenueig.log,cct.council.debt.log ))

cct.council.revenues.log2 <- with(data_rd, rdrobust(y=revenue.PC.D23.log , x=demshare, cluster=cluster, level=90, weight=weight))

cct.council.revenuesown.log2 <- with(data_rd, rdrobust(y=revenueown.PC.D23.log , x=demshare, cluster=cluster, level=90, weight=weight))

 cct.council.taxes.log2 <- with(data_rd, rdrobust(y=Total.Taxes.D23.log, x=demshare, cluster=cluster, level=90, weight=weight))
	
 cct.council.salestaxes.log2 <- with(data_rd, rdrobust(y= salestax.PC.D23.log, x=demshare, cluster=cluster, level=90, weight=weight))
 
  cct.council.salestaxes.select.log2 <- with(data_rd, rdrobust(y=salestax.select.PC.D23.log, x=demshare, cluster=cluster, level=90,weight=weight))
  
 cct.council.propertytaxes.log2 <- with(data_rd, rdrobust(y= propertytax.PC.D23.log, x=demshare, cluster=cluster, level=90, weight=weight))
 
  cct.council.debt.log2 <- with(data_rd, rdrobust(y= debt.PC.D23.log, x=demshare, cluster=cluster, level=90, weight=weight))
  cct.council.charges.log2 <- with(data_rd, rdrobust(y= charges.PC.D23.log, x=demshare, cluster=cluster, weight=weight,level=90))
  cct.council.revenueig.log2 <- with(data_rd, rdrobust(y= revenueig.PC.D23.log, x=demshare, cluster=cluster, weight=weight,level=90))


coefs = data.frame(rbind(cbind(cct.council.revenues.log$coef[1,1],cct.council.revenues.log$ci[3,1],cct.council.revenues.log$ci[3,2],cct.council.revenues.log2$ci[3,1],cct.council.revenues.log2$ci[3,2]),
                         #cbind(cct.council.revenuesown.log$coef[1,1],cct.council.revenuesown.log$ci[3,1],cct.council.revenuesown.log$ci[3,2],cct.council.revenuesown.log2$ci[3,1],cct.council.revenuesown.log2$ci[3,2]),
                         
cbind( cct.council.taxes.log$coef[1,1], cct.council.taxes.log$ci[3,1], cct.council.taxes.log$ci[3,2], cct.council.taxes.log2$ci[3,1], cct.council.taxes.log2$ci[3,2]),
                         cbind(cct.council.salestaxes.log$coef[1,1],cct.council.salestaxes.log$ci[3,1],cct.council.salestaxes.log$ci[3,2],cct.council.salestaxes.log2$ci[3,1],cct.council.salestaxes.log2$ci[3,2]),
                         cbind( cct.council.propertytaxes.log$coef[1,1], cct.council.propertytaxes.log$ci[3,1], cct.council.propertytaxes.log$ci[3,2], cct.council.propertytaxes.log2$ci[3,1], cct.council.propertytaxes.log2$ci[3,2]),
                             cbind( cct.council.charges.log$coef[1,1], cct.council.charges.log$ci[3,1], cct.council.charges.log$ci[3,2], cct.council.charges.log2$ci[3,1], cct.council.charges.log2$ci[3,2]),
                         cbind( cct.council.revenueig.log$coef[1,1], cct.council.revenueig.log$ci[3,1], cct.council.revenueig.log$ci[3,2], cct.council.revenueig.log2$ci[3,1], cct.council.revenueig.log2$ci[3,2]),                                           cbind(cct.council.debt.log$coef[1,1],cct.council.debt.log$ci[3,1],cct.council.debt.log$ci[3,2],cct.council.debt.log2$ci[3,1],cct.council.debt.log2$ci[3,2])

))
names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("Total Revenues",
               "Total Taxes",
               "Sales Taxes",
               "Property Taxes",
               "Charges and Misc. Rev.",
               "Intergov. Rev.",
               "Debt")



####
## Figure 10: The effect of legislative partisanship on changes in logged per capita revenues in the fiscal years two and three years after an election.
####

pdf("figures/figure10_coefplot_revenues_log.pdf",height=5,width=7.5)
ggplot(coefs, aes(c(7:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
     geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                     colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient + Robust CI\n(Ave. in Years 2 and 3 after Election)",breaks=seq(-0.8, .6, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  .3),
                         labels=seq(-0.8, .6, 0.1)) +
      scale_x_continuous("Outcome Variable (logged per capita)",breaks=c(7:1),
                         limits=c(1,7.3),
                         labels=coefs$vars) + 
      coord_flip() + 
        geom_text(nudge_x=.25, size=3) +
    theme_bw()
dev.off()




####
## Table B1: Covariate continuity tests for the County Legislative RD design
## (placebo tests)
##

data_rd<-arrange(data_rd,fips, YearData, district)

data_rd <- group_by(data_rd, fips) %>%
    mutate(
            Total.Expenditure.PC.D6.log = log(Total.Expenditure.PC.L6+.01)-log(Total.Expenditure.PC+.01))
      
data_rd$ demshare.yearlagD<-data_rd$YearData-as.numeric(as.vector(data_rd$ demshare.yearlag))

data.placebocheck<-subset(data_rd, !is.na(Total.Expenditure.PC.D23.log))

data.placebocheck<-subset(data.placebocheck, !is.na(fips))

placebo1 <- with(data.placebocheck[data.placebocheck$demshare.yearlagD<5 & data.placebocheck$demshare.yearlagD>0,], rdrobust(y=demshare.lag1 , x=demshare, weight=weight ))
with(subset(data.placebocheck, demshare>-.1 & demshare<.1), rdplot(y=demshare.lag1 , x=demshare, p=1, nbins=20, y.lim=c(-.2,.2)))

placebo2 <- with(data.placebocheck[data.placebocheck$demshare.yearlagD<5 & data.placebocheck$demshare.yearlagD>0,], rdrobust(y=demseat3 , x=demshare, weight=weight  ))
with(subset(data.placebocheck, demshare>-.1 & demshare<.1), rdplot(y=demseat3 , x=demshare, p=1, nbins=20, y.lim=c(.1,.9)))

placebo3 <- with(data.placebocheck, rdrobust(y=log(Total.Expenditure.PC) , x=demshare ,  weight=weight  ))

placebo4 <- with(data.placebocheck, rdrobust(y=Total.Expenditure.PC.lagged.D1.log , x=demshare, weight=weight ))

placebo5 <- with(data.placebocheck, rdrobust(y=demshare_council2 , x=demshare , weight=weight ))

rd.export(list(placebo1,placebo2,placebo3,placebo4, placebo5))

summary(with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare , weight=weight, cov=demseat2 )))

####
## Figure 7: The effect of legislative partisanship on changes in logged per capita county expenditures in the fiscal 
## years two and three years after an election
####

data.graph<-subset(data_rd, abs(demshare)<.45 & !is.na(demshare) & !is.na(Total.Expenditure.PC.D23.log))
width <- .005
bw.mserd.tri2 <- with(data_rd,
                       rdbwselect(y=Total.Expenditure.PC.D23.log, x=demshare, bwselect="mserd",
                                  kernel='tri', weight=weight))
(mserd.tri2 <- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight)))

data.graph <- mutate(data.graph, bin=cut(demshare, breaks=seq(-1,1, width)))

data.graph2<-group_by(data.graph, bin)
data.graph2<-dplyr::select(data.graph2, bin, Total.Expenditure.PC.D23.log, weight)
data.graph2%>% group_by( bin) %>% 
	summarise_all(funs(sum(. * weight , na.rm = TRUE),  # Weighted Sum
                      weighted.mean(.,w = weight, na.rm = TRUE))) -> agg.df

bin1.df <- data.frame(bin.mean=tapply(data.graph$Total.Expenditure.PC.D23.log,
                         data.graph$bin, mean, na.rm=TRUE, weight=weight),
                     mid=seq(-1 + width/2, 1 - width/2, width),
                     n=tapply(data.graph$Total.Expenditure.PC.D23.log, data.graph$bin, length))
bin1.df$bin<-rownames(bin1.df)
              
bin1.df<-merge(bin1.df, agg.df, by="bin")
cor( bin1.df$bin.mean, bin1.df$Total.Expenditure.PC.D23.log_weighted.mean, use="complete.obs")

library(ggplot2)
pdf("figures/figure7_CouncilRd_Expenditure_logged_weighted.pdf", width=8, height=4)
(ggplot(data.graph)
# + geom_point(aes(x = demshare, y = Total.Expenditure.PC.D23.log, alpha=.1)
 + geom_smooth(data=subset(data.graph, demshare > 0),
               aes(x = demshare, y = Total.Expenditure.PC.D23.log,
                   weight=tri(demshare, bw.mserd.tri2$bws[, "h (left)"])/seats_tot),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, color="black")
 + geom_smooth(data=subset(data.graph, demshare < 0),
               aes(x = demshare, y = Total.Expenditure.PC.D23.log,
                   weight=tri(demshare, bw.mserd.tri2$bws[, "h (left)"])/seats_tot),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, color="black")
 + geom_point(data=bin1.df, aes(x=mid, y=Total.Expenditure.PC.D23.log_weighted.mean, size=n), pch=1)
 + geom_vline(xintercept=0)
 + geom_hline(yintercept=0, lty='dotted')
 + coord_cartesian(ylim = c(-.13, .13))
 + scale_x_continuous(breaks=seq(-1, 1, .025),
                      limits=c(-bw.mserd.tri2$bws[, "h (left)"],
                          bw.mserd.tri2$bws[, "h (left)"]),
                      labels=seq(-100, 100, 2.5))
 + theme_bw()
 +theme(plot.title = element_text(hjust = 0.5))
 + labs(size = '# Obs. in\n0.5% Bin',
        x = "Democratic Margin in County Legislature Election (%)",
        y = "Change in Logged Per Capita Expenditures\n (Ave. of Two Years after Election)")
 + ggtitle("Effect of Partisan Selection on Logged County Expenditures")
 + annotate('text', x = 0.001, y = .1, 0.05, hjust=0, parse=TRUE,
            label = paste('hat(tau)==', round(mserd.tri2$coef['Conventional', ], 3),
                '~(list(', round(mserd.tri2$ci['Robust', 1], 3),
                ',', round(mserd.tri2$ci['Robust', 2], 3), '))'))
 )
dev.off()

####
## Figure 8: The effect of legislative partisanship on per capita county expenditures 1-6 years after an election.
####

## By Time leads
##
cct.expenditures.t1 <- with(data_rd[data_rd$YearData<2011,],
                            rdrobust(y=Total.Expenditure.PC.D1.log ,
                                     x=demshare ))
cct.expenditures.t2 <- with(data_rd[data_rd$YearData<2011,],
                            rdrobust(y=Total.Expenditure.PC.D2.log ,
                                     x=demshare))
cct.expenditures.t3 <- with(data_rd[data_rd$YearData<2011,],
                            rdrobust(y=Total.Expenditure.PC.D3.log ,
                                     x=demshare))
cct.expenditures.t4 <- with(data_rd[data_rd$YearData<2011,],
                            rdrobust(y=Total.Expenditure.PC.D4.log ,
                                     x=demshare))
cct.expenditures.t5 <- with(data_rd[data_rd$YearData<2011,],
                            rdrobust(y=Total.Expenditure.PC.D5.log ,
                                     x=demshare))
cct.expenditures.t6 <- with(data_rd[data_rd$YearData<2011,],
                            rdrobust(y=Total.Expenditure.PC.D6.log ,
                                     x=demshare))


rd<-matrix(nrow=6, ncol=4, data=NA)
rd[,1] <- c(1,2,3,4,5,6)
rd[,2] <- as.numeric(cbind(cct.expenditures.t1$coef['Conventional', ],
                           cct.expenditures.t2$coef['Conventional', ],
                           cct.expenditures.t3$coef['Conventional', ],
                           cct.expenditures.t4$coef['Conventional', ],
                           cct.expenditures.t5$coef['Conventional', ],
                           cct.expenditures.t6$coef['Conventional', ]))

rd[,3] <- as.numeric(cbind(cct.expenditures.t1$ci['Robust', ]["CI Lower"],
                           cct.expenditures.t2$ci['Robust', ]["CI Lower"],
                           cct.expenditures.t3$ci['Robust', ]["CI Lower"],
                           cct.expenditures.t4$ci['Robust', ]["CI Lower"],
                           cct.expenditures.t5$ci['Robust', ]["CI Lower"],
                           cct.expenditures.t6$ci['Robust', ]["CI Lower"]))

rd[,4] <-
    as.numeric(cbind( cct.expenditures.t1$ci['Robust', ]["CI Upper"],
                           cct.expenditures.t2$ci['Robust', ]["CI Upper"],
                           cct.expenditures.t3$ci['Robust', ]["CI Upper"],
                           cct.expenditures.t4$ci['Robust', ]["CI Upper"],
                           cct.expenditures.t5$ci['Robust', ]["CI Upper"],
                           cct.expenditures.t6$ci['Robust', ]["CI Upper"]))
colnames(rd)<-c("lead","est", "ci_lower", "ci_upper")
rd<-as.data.frame(rd)
pdf("figures/figure8_expendituresRD_leads.pdf", width=10, height=4)
(ggplot(subset(rd))
 + geom_pointrange(aes(x = lead, y = est, ymin = `ci_lower`,
                       ymax = `ci_upper`))
 + geom_hline(yintercept = 0, linetype = 'dotted')
 + ylab('Effect on Logged Per Capita\nExpenditures (Robust 95% CI)')
 + xlab("Years after Election")
 + scale_x_continuous(breaks=c(1,2,3,4,5,6),
                      labels=c(1,2,3,4,5,6))
 + theme_bw()
 + ylim(c(-.05,.14))
      + theme(axis.title.x = element_text(size=20),axis.title.y = element_text(size=15),axis.text.y = element_text(size=12),axis.ticks.length=unit(0.5, "cm"))
 )
dev.off()


#####
### Table C1: Robustness of Main Results to Different Modeling Choices.
#####

rd.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster))

rd.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))


rd.bandwidth1.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.01, kernel="uni"))

rd.bandwidth2.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.02, kernel="uni"))

rd.bandwidth3.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.03, kernel="uni"))

rd.bandwidth4.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.04, kernel="uni"))

rd.bandwidth5.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.05, kernel="uni"))

rd.bandwidth1.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.01, weight=weight, kernel="uni"))

rd.bandwidth2.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.02, weight=weight, kernel="uni"))

rd.bandwidth3.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.03, weight=weight, kernel="uni"))

rd.bandwidth4.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.04, weight=weight, kernel="uni"))

rd.bandwidth5.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, h=.05, weight=weight, kernel="uni"))

rd.polynomial1.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=1))

rd.polynomial2.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=2))

rd.polynomial3.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=3))

rd.polynomial4.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=4))

rd.polynomial1.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=1, weight=weight))

rd.polynomial2.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=2, weight=weight))

rd.polynomial3.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=3, weight=weight))

rd.polynomial4.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, p=4, weight=weight))

## manually change confidence intervals for hard-wired bandwidths to non-optimal p-values and confidence intervals
### Manually change labels
rd.export(list(rd.weighted,rd.bandwidth1.weighted, rd.bandwidth2.weighted, rd.bandwidth3.weighted, rd.bandwidth4.weighted, rd.bandwidth5.weighted, rd.polynomial2.weighted,rd.polynomial3.weighted,rd.polynomial4.weighted))

rd.export(list(rd.nonweighted,rd.bandwidth1.nonweighted, rd.bandwidth2.nonweighted, rd.bandwidth3.nonweighted, rd.bandwidth4.nonweighted, rd.bandwidth5.nonweighted, rd.polynomial2.nonweighted,rd.polynomial3.nonweighted,rd.polynomial4.nonweighted))

 ### Randomization inference
tmp <- with(data_rd,rdrandinf(Y= Total.Expenditure.PC.D23.log , R=  demshare,wl=-.01,wr=.01, rep=10000))

#####
## Analysis of CF Scores for Figure 5: The effect of legislative partisanship on legislator CF-Score.
#####

council_all<-read.dta(file="council_all.dta")
county_elections_cfscores<-read.csv(file="county_elections_cfscores.csv")
county_elections_cfscores<-subset(county_elections_cfscores, variable=="cand1")

## multimember districts
county_elections_cfscores_mm<-read.csv(file="county_elections_cfscores_mm.csv")
council_all_mm<-dplyr::select(council_all, fips, year, district, n_winners)
county_elections_cfscores_mm<-merge(county_elections_cfscores_mm, council_all_mm, by=c("fips", "year", "district"))

county_elections_cfscores_mm2<-subset(county_elections_cfscores_mm, n_winners==2)
county_elections_cfscores_mm2<-subset(county_elections_cfscores_mm2, variable=="cand1" | variable=="cand2")

county_elections_cfscores_mm3<-subset(county_elections_cfscores_mm, n_winners==3)
county_elections_cfscores_mm3<-subset(county_elections_cfscores_mm3, variable=="cand1" | variable=="cand2" | variable=="cand3")

county_elections_cfscores_mm4<-subset(county_elections_cfscores_mm, n_winners==4)
county_elections_cfscores_mm4<-subset(county_elections_cfscores_mm4, variable=="cand1" | variable=="cand2" | variable=="cand3"| variable=="cand4" )

county_elections_cfscores_mm5<-subset(county_elections_cfscores_mm, n_winners==5)
county_elections_cfscores_mm5<-subset(county_elections_cfscores_mm5, variable=="cand1" | variable=="cand2" | variable=="cand3"| variable=="cand4" | variable=="cand5")

county_elections_cfscores_mm<-rbind(county_elections_cfscores_mm2, county_elections_cfscores_mm3, county_elections_cfscores_mm4, county_elections_cfscores_mm5)

council_all<-read.dta(file="council_all.dta")
council_all<-merge(council_all, county_elections_cfscores, by=c("fips", "year", "district"))

council_all<-subset(council_all, n_winners==1)
council_all<-merge(council_all, council_all_mm, all.x=T, all.y=T)

council_means_parties<-subset(council_all, variable=="cand1")
council_means_parties<-subset(council_means_parties, !is.na(cfscore))
council_means_parties <- group_by(council_means_parties, fips,winner_pid) %>%
	summarise(
    median_cfscore_party = mean(cfscore, na.rm=T),  
    length_cfscore_party = length(cfscore)  
    )
council_means_parties<-subset(council_means_parties, length_cfscore_party>2)
## Council composition data. Has total number of seats in each council-year.
councilcomp_export<-read.csv("councilcomp_export.csv") 
councilcomp_export<-dplyr::select(councilcomp_export, -population_est)
## this file has the number of seats in counties that are missing from our panel
council_seats<-read.csv("counties_missing_seats_total.csv")
# colnames(council_seats) <- c("X","fips")

councilcomp_export<-merge(councilcomp_export, council_seats, by.x=c( "fips"), by.y=c("fips"), all.x=T)

council_all<-merge(council_all, councilcomp_export, by.x=c( "fips","year"), by.y=c("fips","year"), all.x=T, all.y=T)
council_all$ year<-as.numeric(as.vector(council_all$ year))
council_all$partisan_diff<-abs(council_all$seats_dem-council_all$seats_rep)

council_all<-arrange(council_all, fips, year, district)
council_all <- group_by(council_all, fips) %>%
    mutate(
    partisan_diff.L1 = lead.new(partisan_diff, 1, along_with = year),
      seats_dem.L1 = lead.new(seats_dem, 1, along_with = year),
        seats_rep.L1 = lead.new(seats_rep, 1, along_with = year),
          seats_tot.L1 = lead.new(seats_tot, 1, along_with = year))


council_all$winner_dem<-NA
council_all$winner_dem[council_all$demshare>0]<-1
council_all$winner_dem[council_all$demshare<0]<-0

council_all$seats_dem.L1.counterfactual<-council_all$seats_dem.L1
council_all$seats_rep.L1.counterfactual<-council_all$seats_rep.L1

council_all$seats_dem.L1.counterfactual[council_all$winner_dem==1& !is.na(council_all$winner_dem)]<- council_all$seats_dem.L1.counterfactual[council_all$winner_dem==1 & !is.na(council_all$winner_dem)]-1

council_all$seats_rep.L1.counterfactual[council_all$winner_dem==0& !is.na(council_all$winner_dem)]<-council_all$seats_rep.L1.counterfactual[council_all$winner_dem==0& !is.na(council_all$winner_dem)]-1

council_all$counterfactual_tied<-0

council_all$counterfactual_partisandiff<-abs(council_all$seats_dem.L1.counterfactual-council_all$seats_rep.L1.counterfactual)
table(council_all$counterfactual_partisandiff, council_all$partisan_diff)


### re-weight the data
council_all$mean_seats_tot<-mean(council_all$seats_tot, na.rm=T)
council_all$weight<-council_all$mean_seats_tot/council_all$seats_tot
council_all<-group_by(council_all, fips, year) %>%
	mutate(sum_weights=sum(weight))

years<-c(1992,1994, 1996, 1998, 2000, 2002, 2004 , 2006, 2008, 2010, 2012)
council_median<-NA
for (i in 1:length(years)){
	sub<-subset(council_all, year==years[i])
	sub<-dplyr::select(sub, fips, district, year, n_winners, cand1, winner_pid)
	sub2<-subset(council_all, year==years[i]-2)
	sub2<-dplyr::select(sub2, fips, district, year, n_winners, cand1, winner_pid)
	sub3<-rbind(sub, sub2)
	cfscores_unique<-group_by(county_elections_cfscores, fips,district, value, variable) %>%
		summarise(cfscore=first(cfscore))
	sub4<-merge(sub3, cfscores_unique,by.x=c("fips",  "district", "cand1"), by.y=c("fips",  "district", "value"))
	
	sub4$duplicate<-duplicated(paste(sub4$fips, sub4$district, sub4$cand1))
	sub4<-subset(sub4, duplicate==F)
	sub4<-merge(sub4, council_means_parties, by=c("fips", "winner_pid"), all.x=T)
	sub4$cfscore[is.na(sub4$cfscore)]<-sub4$median_cfscore_party[is.na(sub4$cfscore)]
	#sub4<-subset(sub4, !is.na(cfscore))
	sub5 <- group_by(sub4, fips) %>%
	summarise(
    median_cfscore = median(cfscore, na.rm=T),
    length_cfscores = length(cfscore)
   )
   sub5$year<-years[i]
   council_median<-rbind(council_median, sub5)
}


council_all<-merge(council_all, council_median, by=c("year", "fips"), all.x=T)
council_all$cfscores_availability<-council_all$length_cfscores/council_all$seats_tot

council_all<-arrange(council_all, fips, year, district)


rd<- with(council_all, rdrobust(y=cfscore , x=demshare, weight=weight ))

####
## Figure 5: The effect of legislative partisanship on legislator CF-Score.
####
council_graph<-council_all
tri <- function (x, h, c=0) pmax(0, 1 - abs((x - c) / h))
width <- .005
bw.mserd.tri2 <- with(council_graph,
                       rdbwselect(y=cfscore, x=demshare, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri2 <- with(council_graph, rdrobust(y=cfscore, x=demshare, weight=weight)))
council_graph<-subset(council_graph, !is.na(cfscore))


council_graph <- mutate(council_graph, bin=cut(demshare, breaks=seq(-1,1, width)))
bin1.df <- data.frame(bin.mean=tapply(council_graph$cfscore,
                         council_graph$bin, mean, na.rm=TRUE),
                     mid=seq(-1 + width/2, 1 - width/2, width),
                     n=tapply(council_graph$cfscore, council_graph$bin, length))
bin1.df <-subset(bin1.df , mid>-.15 & mid<.15)

pdf("Figures/figure5_cfscores.pdf", width=7.5, height=4)
(ggplot(council_graph)
# + geom_point(aes(x = demshare, y = cfscore), alpha=.1)
 + geom_point(data=council_graph, aes(x=demshare, y=cfscore), pch=16, color="grey", size=1)
 + geom_smooth(data=subset(council_graph, demshare > 0),
               aes(x = demshare, y = cfscore,
                   weight=tri(demshare, bw.mserd.tri2$bws[, "h (left)"])),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, color="black")
 + geom_smooth(data=subset(council_graph, demshare < 0),
               aes(x = demshare, y = cfscore,
                   weight=tri(demshare, bw.mserd.tri2$bws[, "h (left)"])),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5, color="black")
 + geom_point(data=bin1.df, aes(x=mid, y=bin.mean, size=n), pch=1)
 + geom_vline(xintercept=0)
 + geom_hline(yintercept=0, lty='dotted')
 + coord_cartesian(ylim = c(-1.5, 1.5))
 + scale_x_continuous(breaks=seq(-1, 1, .05),
                      limits=c(-bw.mserd.tri2$bws[, "h (left)"],
                          bw.mserd.tri2$bws[, "h (left)"]),
                      labels=seq(-100, 100, 5))
 + theme_bw()
 +theme(plot.title = element_text(hjust = 0.5))
 + labs(size = '# Obs. in\n0.5% Bin',
        x = "Democratic Margin in County Legislature Election (%)",
        y = "Legislator CF-Score")
 + ggtitle("Effect of Partisan Selection on County Legislator CF-Scores")
 + annotate('text', x = 0.001, y = 1, 0.05, hjust=0, parse=TRUE,
            label = paste('hat(tau)==', round(mserd.tri2$coef['Conventional', ], 3),
                '~(list(', round(mserd.tri2$ci['Robust', 1], 3),
                ',', round(mserd.tri2$ci['Robust', 2], 3), '))'))
 )
dev.off()


##### 
##### Pivotal legislature for majority control analysis
##### 

data_rd$winner_dem<-NA
data_rd$winner_dem[data_rd$demshare>0]<-1
data_rd$winner_dem[data_rd$demshare<0]<-0

data_rd$seats_dem.L1.counterfactual<-data_rd$seats_dem.L1
data_rd$seats_rep.L1.counterfactual<-data_rd$seats_rep.L1

data_rd$seats_dem.L1.counterfactual[data_rd$winner_dem==1 & !is.na(data_rd$winner_dem)]<- data_rd$seats_dem.L1.counterfactual[data_rd$winner_dem==1 & !is.na(data_rd$winner_dem)]-1

data_rd$seats_rep.L1.counterfactual[data_rd$winner_dem==0& !is.na(data_rd$winner_dem)]<-data_rd$seats_rep.L1.counterfactual[data_rd$winner_dem==0& !is.na(data_rd$winner_dem)]-1

data_rd$counterfactual_tied<-0

data_rd$counterfactual_demshare<-data_rd$seats_dem.L1.counterfactual/(data_rd$seats_dem.L1.counterfactual+data_rd$seats_rep.L1.counterfactual)

data_rd$counterfactual_partisandiff<-abs(data_rd$seats_dem.L1.counterfactual-data_rd$seats_rep.L1.counterfactual)
table(data_rd$counterfactual_partisandiff, data_rd$partisan_diff)


#####
### Moderators of results - main text
#####

writeLines(prettyNum(nrow(data_rd),big.mark=","),"Tables/n_elecs.tex")

rd.nonweighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster))

rd.weighted <-  with(data_rd,
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

####
## Figure 11: The effect of legislative partisanship on changes in logged per capita expenditures by the partisan majority’s margin.
####

cct.council.notclose <- with(data_rd[(data_rd$counterfactual_partisandiff>2),], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))
cct.council.close <- with(data_rd[(data_rd$counterfactual_partisandiff<3),], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))
rd.export(list(cct.council.notclose,cct.council.close))


close_rd <- (with(data_rd[(data_rd$counterfactual_partisandiff<3) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight)))
close_rd_90 <- (with(data_rd[(data_rd$counterfactual_partisandiff<3) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight,level = 90)))

extreme_rd <- (with(data_rd[(data_rd$counterfactual_partisandiff>2) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight)))
extreme_rd_90 <- (with(data_rd[(data_rd$counterfactual_partisandiff>2) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight,level=90)))

coefs = data.frame(rbind(cbind(close_rd$coef[1,1],close_rd$ci[3,1],close_rd$ci[3,2],close_rd_90$ci[3,1],close_rd_90$ci[3,2]),
                         cbind(extreme_rd$coef[1,1],extreme_rd$ci[3,1],extreme_rd$ci[3,2],extreme_rd_90$ci[3,1],extreme_rd_90$ci[3,2])
))

names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("<3 Seat Margin",
               "3+ Seat Margin")

rd.export(list(close_rd,extreme_rd))

pdf("Figures/figure11_coefplot_exp_log_majoritymargin.pdf",height=2,width=7)
ggplot(coefs, aes(c(2:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
      geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                    colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient on Logged Per Capita Spending + Robust CI\n(Ave. in Years 2 and 3 after Election)",
                         breaks=seq(-0.4, .8, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  max(coefs$ci95.hi)),
                         labels=seq(-0.4, .8, 0.1)) +
      scale_x_continuous("Majority Party Margin",
                         breaks=c(2:1),
                         limits=c(0.7,2.3),
                         labels=coefs$vars,position = "bottom") + 
      coord_flip() + 
      geom_text(nudge_x=.2, size=3) +
      theme_bw()
dev.off()




####
## Figure 12: The effect of legislative partisanship on changes in logged per capita expenditures by the size of the county legislature.
####

rd.weighted.3 <-  with(data_rd[data_rd$seats_tot==3,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.4.10 <-  with(data_rd[data_rd$seats_tot>3 & data_rd$seats_tot<11,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.11 <-  with(data_rd[data_rd$seats_tot>10,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.export(list(rd.weighted.3,rd.weighted.4.10,rd.weighted.11))

cct.council.exp.log0 <- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

cct.council.exp.log1 <- with(data_rd[data_rd$seats_tot==3,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))
cct.council.exp.log2 <- with(data_rd[data_rd$seats_tot>3 & data_rd$seats_tot<11,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))
cct.council.exp.log3 <- with(data_rd[data_rd$seats_tot>10  ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))
#cct.council.exp.log4 <- with(data_rd[data_rd$seats_tot>19,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight))
rd.export(list(cct.council.exp.log0,cct.council.exp.log1,cct.council.exp.log2,cct.council.exp.log3))

cct.council.exp.log02 <- with(data_rd, rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90))
cct.council.exp.log12 <- with(data_rd[data_rd$seats_tot==3,], rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight, level=90))
cct.council.exp.log22 <- with(data_rd[data_rd$seats_tot>3 & data_rd$seats_tot<11,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90))
cct.council.exp.log32 <- with(data_rd[data_rd$seats_tot>10  ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90))
#cct.council.exp.log42 <- with(data_rd[data_rd$seats_tot>19,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight, level=90))

coefs = data.frame(rbind(cbind(cct.council.exp.log0$coef[1,1],cct.council.exp.log0$ci[3,1],cct.council.exp.log0$ci[3,2],cct.council.exp.log02$ci[3,1],cct.council.exp.log02$ci[3,2]),
                         cbind(cct.council.exp.log1$coef[1,1],cct.council.exp.log1$ci[3,1],cct.council.exp.log1$ci[3,2],cct.council.exp.log12$ci[3,1],cct.council.exp.log12$ci[3,2]),
                         
cbind( cct.council.exp.log2$coef[1,1], cct.council.exp.log2$ci[3,1], cct.council.exp.log2$ci[3,2], cct.council.exp.log22$ci[3,1], cct.council.exp.log22$ci[3,2]),
                         cbind(cct.council.exp.log3$coef[1,1],cct.council.exp.log3$ci[3,1],cct.council.exp.log3$ci[3,2],cct.council.exp.log32$ci[3,1],cct.council.exp.log32$ci[3,2])
                         
))

names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("Average County",
               "Counties with 3 Legislators",
               "Counties with 4-10 Legislators",
               "Counties with 11+ Legislators")

pdf("figures/figure12_coefplot_size_log.pdf",height=4,width=7.5)
ggplot(coefs, aes(c(4:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
     geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                     colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient + Robust CI\n(Ave. of County Expenditures in Years 2 and 3 after Election)",breaks=seq(-0.8, .6, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  .4),
                         labels=seq(-0.8, .6, 0.1)) +
      scale_x_continuous("",breaks=c(4:1),
                         limits=c(1,4.3),
                         labels=coefs$vars) + 
      coord_flip() + 
        geom_text(nudge_x=.15, size=4) +
    theme_bw()
dev.off()


####
## Figure 13: The effect of legislative partisanship on changes in logged per capita expenditures by the county’s form of government.
####
rd.weighted.fog1 <-  with(data_rd[data_rd$fog_combined=="1",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.fog2 <-  with(data_rd[data_rd$fog_combined=="2",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.fog3 <-  with(data_rd[data_rd$fog_combined=="3",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.export(list(rd.weighted.fog1,rd.weighted.fog2,rd.weighted.fog3))


commission_rd <- (with(data_rd[(data_rd$fog_combined==1) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight)))
commission_rd_90 <- (with(data_rd[(data_rd$fog_combined==1) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight,level = 90)))

council_rd <- (with(data_rd[(data_rd$fog_combined==2) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight)))
council_rd_90 <- (with(data_rd[(data_rd$fog_combined==2) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight,level = 90)))

exec_rd <- (with(data_rd[(data_rd$fog_combined==3) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight)))
exec_rd_90 <- (with(data_rd[(data_rd$fog_combined==3) ,], rdrobust(y=Total.Expenditure.PC.D23.log , x=demshare, cluster=cluster, weight=weight,level = 90)))


coefs = data.frame(rbind(cbind(commission_rd$coef[1,1],commission_rd$ci[3,1],commission_rd$ci[3,2],commission_rd_90$ci[3,1],commission_rd_90$ci[3,2]),
                         cbind(council_rd$coef[1,1],council_rd$ci[3,1],council_rd$ci[3,2],council_rd_90$ci[3,1],council_rd_90$ci[3,2]),
                         cbind(exec_rd$coef[1,1],exec_rd$ci[3,1],exec_rd$ci[3,2],exec_rd_90$ci[3,1],exec_rd_90$ci[3,2])
))

names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c("Commission",
               "Council-manager",
               "Council + Elected Executive")

rd.export(list(commission_rd,council_rd,exec_rd))

pdf("Figures/figure13_coefplot_exp_log_fog.pdf",height=3,width=7)
ggplot(coefs, aes(c(3:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
      geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                    colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient on Logged Per Capita Spending + Robust CI\n(Ave. in Years 2 and 3 after Election)",
                         breaks=seq(-0.4, .8, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  max(coefs$ci95.hi)),
                         labels=seq(-0.4, .8, 0.1)) +
      scale_x_continuous("Form of Government",
                         breaks=c(3:1),
                         limits=c(0.8,3.2),
                         labels=coefs$vars,position = "bottom") + 
      coord_flip() + 
      geom_text(nudge_x=.2, size=3) +
      theme_bw()
dev.off()


####
## Figure 14: The effect of legislative partisanship on changes in logged per capita expenditures by the percent of revenue from intergovernmental sources. 
####
## variation in ig rev by state
igrev<-group_by(data.use[data.use$population_est>150000,], STATEAB)  %>%
	summarize(igrev.perc.state=mean(igrev.PC)/mean(revenue.PC))
data_rd<-merge(data_rd, igrev, by="STATEAB", all.x=T)
ig1<-with(data_rd[data_rd$igrev.perc.state>.35,],
					rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
ig2<-with(data_rd[data_rd$igrev.perc.state<.35,],
					rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.export(list(ig1, ig2))
rd.export(list(ig1, ig2))



ighi_rd <- with(data_rd[data_rd$igrev.perc.state>.35,],
          rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
ighi_rd_90 <- with(data_rd[data_rd$igrev.perc.state>.35,],
                rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight,level = 90))

iglo_rd <- with(data_rd[data_rd$igrev.perc.state<.35,],
          rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
iglo_rd_90 <- with(data_rd[data_rd$igrev.perc.state<.35,],
                rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight,level = 90))


coefs = data.frame(rbind(cbind(ighi_rd$coef[1,1],ighi_rd$ci[3,1],ighi_rd$ci[3,2],ighi_rd_90$ci[3,1],ighi_rd_90$ci[3,2]),
                         cbind(iglo_rd$coef[1,1],iglo_rd$ci[3,1],iglo_rd$ci[3,2],iglo_rd_90$ci[3,1],iglo_rd_90$ci[3,2])
))

names(coefs) = c("coef","ci95.lo","ci95.hi","ci90.lo","ci90.hi") 
coefs$vars = c(">35% of Revenue",
               "<35% of Revenue")

rd.export(list(ighi_rd,iglo_rd))

pdf("Figures/figure14_coefplot_exp_log_igrev.pdf",height=2,width=7)
ggplot(coefs, aes(c(2:1), coef, label=round(coef,2)))+
      geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
      geom_errorbar(aes(ymin=ci95.lo, ymax=ci95.hi), 
                    colour="black", width=0, size=.5) +
      geom_errorbar(aes(ymin=ci90.lo, ymax=ci90.hi), 
                    colour="black", width=0, size=1) +
      geom_point(size=4, pch=21, fill="black") +
      scale_y_continuous("RD Coefficient on Logged Per Capita Spending + Robust CI\n(Ave. in Years 2 and 3 after Election)",
                         breaks=seq(-0.4, .8, 0.1),
                         limits=c(min(coefs$ci95.lo),
                                  max(coefs$ci95.hi)),
                         labels=seq(-0.4, .8, 0.1)) +
      scale_x_continuous("% of Revenue from\nIntergovernmental\nSources",
                         breaks=c(2:1),
                         limits=c(0.7,2.3),
                         labels=coefs$vars,position = "bottom") + 
      coord_flip() + 
      geom_text(nudge_x=.2, size=3) +
      theme_bw()
dev.off()


#####
### Moderators of results - appendix
#####

####
## Table E2: Heterogeneity in effect of legislator partisanship across regions
####

rd.weighted.nonsouth <-  with(data_rd[ data_rd$Census.Region!=3,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.south <-  with(data_rd[ data_rd$Census.Region==3,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

data_rd$Census.Region2<-NA
data_rd$Census.Region2[data_rd$Census.Region==3]<-"south"
data_rd$Census.Region2[data_rd$Census.Region==1]<-"east"
data_rd$Census.Region2[data_rd$Census.Region==2]<-"midwest"
data_rd$Census.Region2[data_rd$Census.Region==4]<-"west"
rd.weighted.east <-  with(data_rd[data_rd$Census.Region2=="east",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.midwest <-  with(data_rd[data_rd$Census.Region2=="midwest",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.south <-  with(data_rd[data_rd$Census.Region2=="south",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.west <-  with(data_rd[data_rd$Census.Region2=="west",],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.export(list(rd.weighted.east,rd.weighted.midwest,rd.weighted.south,rd.weighted.west))


####
## Table E3: Heterogeneity in effect of legislator partisanship across time
####
rd.weighted.1990s <-  with(data_rd[data_rd$YearData<2003,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.weighted.2000s <-  with(data_rd[data_rd$YearData>2002,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))
rd.export(list(rd.weighted.1990s,rd.weighted.2000s))

####
## Table E4: Heterogeneity in effect of legislator partisanship across more urban and rural counties
####
County_Rural_Lookup<-read.csv(file="County_Rural_Lookup.csv")
data_rd<-merge(data_rd, County_Rural_Lookup, by="fips", all.x=T)

rd.weighted.rural <-  with(data_rd[data_rd$pop_2010_percent_rural>15,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.weighted.urban <-  with(data_rd[data_rd$pop_2010_percent_rural<15.01,],
 	rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.export(list(rd.weighted.rural,rd.weighted.urban))

####
## Table E6: Heterogeneity in effect of legislator partisanship by county size
####
rd.weighted.75k <-  with(data_rd_100kthreshold[data_rd_100kthreshold$population_est>75000 & data_rd_100kthreshold$population_est<150000,],
												 rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.weighted.250k <-  with(data_rd[data_rd$population_est>150000& data_rd$population_est<250000,],
													rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.weighted.400k <-  with(data_rd[data_rd$population_est<400000 & data_rd$population_est>250000,], rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.weighted.750k <-  with(data_rd[data_rd$population_est<800000 & data_rd$population_est>400000,],rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.weighted.big <-  with(data_rd[data_rd$population_est>800000 ,],rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster, weight=weight))

rd.export(list(rd.weighted, rd.weighted.75k,rd.weighted.250k, rd.weighted.400k,rd.weighted.750k,rd.weighted.big))

####
## Table E7: Heterogeneity in effect of legislator partisanship based on partisan control of state government
####
load(file="party_control_data_150710.RData")
data_rd<-merge(data_rd, data2, by.x=c("YearData", "STATEAB"), by.y=c("year", "abb"), all.x=T)

data_rd$state_gov<-NA
data_rd$state_gov<-data_rd$sen_dem_control + data_rd$hs_dem_control+ data_rd$gov_party
data_rd$state_gov<-data_rd$state_gov/3
dem_control<-(with(data_rd[data_rd$state_gov==1,],
									 rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster,  weight=weight)))
rep_control<-(with(data_rd[data_rd$state_gov==0,],
									 rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster,  weight=weight)))
split_control<-(with(data_rd[data_rd$state_gov>0 & data_rd$state_gov<1,],
										 rdrobust(y=Total.Expenditure.PC.D23.log, x=demshare, cluster=cluster,  weight=weight)))

summary(dem_control)
summary(rep_control)
summary(split_control)
rd.export(list(dem_control,rep_control,split_control))

## variation in council size by region, institution
collapsed<-group_by(data_rd,fips,Census.Region,STATEAB) %>%
	summarize(seats_tot=mean(seats_tot)	,
		fog_combined=max(fog_combined, na.rm=T),
		pop_2010_percent_rural=mean(pop_2010_percent_rural, na.rm=T)
	)
collapsed$Census.Region2<-NA
collapsed$seats_tot_bin<-NA
collapsed$Census.Region2[collapsed$Census.Region==3]<-"south"
collapsed$Census.Region2[collapsed$Census.Region==1]<-"east"
collapsed$Census.Region2[collapsed$Census.Region==2]<-"midwest"
collapsed$Census.Region2[collapsed$Census.Region==4]<-"west"
collapsed$seats_tot_bin[collapsed$seats_tot==3]<-"small"
collapsed$seats_tot_bin[collapsed$seats_tot>3 & collapsed$seats_tot<10]<-"medium"
collapsed$seats_tot_bin[collapsed$seats_tot>9]<-"large"

collapsed$seats_tot_bin <- factor(collapsed$seats_tot_bin, levels = c("small", "medium", "large"))


collapsed$urban<-NA
collapsed$urban[collapsed$pop_2010_percent_rural>15]<-"rural"
collapsed$urban[collapsed$pop_2010_percent_rural<15.01]<-"urban"

## Table E5: Proportion of urban and rural counties with small, medium, and large legislatures
xtable(prop.table(table(collapsed$urban,collapsed$seats_tot_bin),1))

##Table E1: Proportion of counties in each region with small, medium, and large legislatures
xtable(prop.table(table(collapsed$Census.Region2,collapsed$seats_tot_bin),1))

##Extra table: Proportion of counties in each FOG with small, medium, and large legislatures
xtable(prop.table(table(collapsed$fog_combined,collapsed$seats_tot_bin),1))
